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ANALYSIS OF MIXED INTERIOR PENALTY DISCONTINUOUS 
GALERKIN METHODS FOR THE CAHN-HILLIARD EQUATION 
AND THE HELE-SHAW FLOW 

XIAOBING FENG*, YUKUN Lit, and YULONG XING* 

Abstract. This paper proposes and analyzes two fully discrete mixed interior penalty discontin¬ 
uous Galerkin (DG) methods for the fourth order nonlinear Cahn-Hilliard equation. Both methods 
use the backward Euler method for time discretization and interior penalty discontinuous Galerkin 
methods for spatial discretization. They differ from each other on how the nonlinear term is treated, 
one of them is based on fully implicit time-stepping and the other uses the energy-splitting time¬ 
stepping. The primary goal of the paper is to prove the convergence of the numerical interfaces of the 
DG methods to the interface of the Hele-Shaw flow. This is achieved by establishing error estimates 
that depend on only in some low polynomial orders, instead of exponential orders. Similar to 
mi, the crux is to prove a discrete spectrum estimate in the discontinuous Galerkin finite element 
space. However, the validity of such a result is not obvious because the DG space is not a subspace of 
the (energy) space H^{Q) and it is larger than the finite element space. This difficult is overcome by 
a delicate perturbation argument which relies on the discrete spectrum estimate in the finite element 
space proved in m Numerical experiment results are also presented to gauge the theoretical results 
and the performance of the proposed fully discrete mixed DG methods. 

Key words. Cahn-Hilliard equation, Hele-Shaw problem, phase transition, discontinuous Galerkin 
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1. Introduction. This paper concerns with mixed interior penalty discontinu- 
ous Galerkin (MIP-DG) approximations of the following Cahn-Hilliard problem: 


( 1 . 1 ) 

( 1 . 2 ) 

(1.3) 

(1.4) 


ut — Aw 

—eAu + -f{u) 
e 

du dw 
dn dn 
u 


= 0 in := H X (0,T), 

= w in 

= 0 on d^T '■= dfl X (0, T), 
= uo in n X {t = 0}. 


Here H C R^ (d = 2,3) is a bounded domain, and f{u) = F'{u), F{u) is a nonconvex 
potential density function which takes its global minimum zero at u = ±1. In this 
paper, we only consider the following quartic potential density function: 

(1.5) F{u) = \{u^-ir. 

After eliminating the intermediate variable w (called the chemical potential), the 
above system reduces into a fourth order nonlinear PDE for u, which is known as the 
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Cahn-Hilliard equation in the literature. This equation was originally introduced by 
John W. Cahn and John E. Hilliard in to describe the process of phase separa¬ 
tion, by which the two components of a binary fluid spontaneously separate and form 
domains pure in each component. Here u and 1 — u denote respectively the concen¬ 
trations of the two fluids, with u = ±1 indicating domains of the two components. 
We note that the equation ([iTf-Ol differs from the original Cahn-Hilliard equation 
in the scaling of the time, and t here corresponds to f in the original formulation, e, 
which is positively small, is called the interaction length. 

Besides its important role in materials phase transition, the Cahn-Hilliard equa¬ 
tion has been extensively studied due to its close relation with the Hele-Shaw prob¬ 
lem. It was first formally proved by Pego [19] that the chemical potential w := 
—eAu + ^f{u) tends to a limit which satisfies the following free boundary problem 
known as the Hele-Shaw problem: 


( 1 . 6 ) 

(1.7) 

( 1 . 8 ) 

(1.9) 


Aw = 0 
dw 

w = an 




dw 


-dn 


inn\Tu te [o,r], 

on t e [0, T], 
on Pt, t e [o,r], 
on Pt, t e [o,r]. 


as e \ 0, provided that the Hele-Shaw problem has a unique classical solution. Here 



K and V represent the mean curvature and the normal velocity of the interface Pj. A 
rigorous justification that u —)■ ±1 in the interior or exterior of P* for all t G [0,T] as 
e \ 0 was given by Stoth for the radially symmetric case, and by Alikakos, Bates 
and Chen |2| for the general case. In addition, Chen [7] established the convergence 
of the weak solution of the Cahn-Hilliard problem to a weak (or varifold) solution of 
the Hele-Shaw problem. 

Moreover, the Cahn-Hilliard equation (together with the Allen-Cahn equation) 
has become a fundamental equation as well as a building block in the phase field 
methodology (or the diffuse interface methodology) for moving interface and free 
boundary problems arising from various applications such as fluid dynamics, mate¬ 
rials science, image processing and biology (cf. [501 [H] and the references therein). 
The diffuse interface approach provides a convenient mathematical formalism for nu¬ 
merically approximating the moving interface problems because explicitly tracking 
the interface is not needed in the diffuse interface formulation. The main advantage 
of the diffuse interface method is its ability to handle with ease singularities of the 
interfaces. Like many singular perturbation problems, the main computational issue 
is to resolve the (small) scale introduced by the parameter e in the equation. Com¬ 
putationally, the problem could become intractable, especially in three-dimensional 
cases if uniform meshes are used. This difficulty is often overcome by exploiting the 
predictable (at least for small e) PDE solution profile and by using adaptive mesh 
techniques (cf. [T6| and the references therein), so fine meshes are only used in the 
diffuse interface region. 

Numerical approximations of the Cahn-Hilliard equation have been extensively 
carried out in the past thirty years (cf. [HI [HI [Hj and the references therein). On 
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the other hand, the majority of these works were done for a fixed parameter e. The 
error bounds, which are obtained using the standard Gronwall inequality technique, 
show an exponential dependence on 1/e. Such an estimate is clearly not useful for 
small e, in particular, in addressing the issue whether the computed numerical inter¬ 
faces converge to the original sharp interface of the Hele-Shaw problem. Better and 
practical error bounds should only depend on 1/e in some (low) polynomial orders 
because they can be used to provide an answer to the above convergence question, 
which in fact is the best result (in terms of e) one can expect. The first such poly¬ 
nomial order in 1/e a priori estimate was obtained in [15] for mixed finite element 
approximations of the Cahn-Hilliard problem (1.11-( 1.51. In addition, polynomial or¬ 


der in 1/e a posteriori error estimates were obtained in |16| for the same mixed finite 
element methods. One of the key ideas employed in all these works is to use a non¬ 
standard error estimate technique which is based on establishing a discrete spectrum 
estimate (using its continuous counterpart) for the linearized Cahn-Hilliard operator. 
An immediate corollary of the polynomial order in 1/e a priori and a posteriori error 
estimates is the convergence of the numerical interfaces of the underlying mixed finite 
element approximations to the Hele-Shaw flow before the onset of singularities of the 
Hele-Shaw flow as e and mesh sizes h and k all tend to zero. 

The objectives of this paper are twofold: Firstly, we develop some MIP-DG meth¬ 
ods and to establish polynomial order in 1/e a priori error bounds, as well as to prove 
convergence of numerical interfaces for the MIP-DG methods. This goal is motivated 
by the advantages of DG methods in regard to designing adaptive mesh methods and 
algorithms, which is an indispensable strategy with the diffuse interface methodology. 
Secondly, we use the Cahn-Hilliard equation as another prototypical model problem 
[18| to develop new analysis techniques for analyzing convergence of numerical inter¬ 
faces to the underlying sharp interface for DG (and nonconforming finite element) 
discretizations of phase field models. To the best of our knowledge, no such con¬ 
vergence result and analysis technique is available in the literature for fourth order 
PDFs. The main obstacle for improving the finite element techniques of m is that 
the DG (and nonconforming finite element) spaces are not subspaces of As a 

result, whether the needed discrete spectrum estimate holds becomes a key question 
to answer. 

This paper consists of four additional sections. In section we first collect some 
a priori error estimates for problem (I.1|-(I.5|, which show the explicit dependence 
on the parameter e. We then cite two important technical lemmas to be used in 
the later sections. One of the lemma states the spectral estimate for the linearized 
Gahn-Hilliard operator. In section we propose two fully discrete MIP-DG schemes 
for problem (l.I)-(1.51, they differ only in their treatment of the nonlinear term. The 


first main result of this section is to establish a discrete spectrum estimate in the DG 
space, which mimics the spectral estimates for the differential operator and its finite 
element counterpart. The second main result of this section is to derive optimal error 
bounds which depends on 1/e only in low polynomial orders for both fully discrete 
MIP-DG methods. In section]^ using the refined error estimates of section]^ we prove 
the convergence of the numerical interfaces of the fully discrete MIP-DG methods to 
the interface of the Hele-Shaw flow before the onset of the singularities as e, h and k 
all tend to zero. Finally, in sectionj^we provide some numerical experiments to gauge 
the performance of the proposed fully discrete MIP-DG methods. 

2. Preliminaries. In this section, we shall collect some known results about 
problem (1.1)-(1.5) from mung, which will be used in sections and Some 
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general assumptions on the initial condition, as well as some energy estimates based on 
these assumptions, will be cited. Standard function and space notations are adopted 
in this paper mm- We use (•,•) and || • 11^2 to denote the standard inner product 
and norm on L^(r2). Throughout this paper, C denotes a generic positive constant 
independent of e, space and time step sizes h and k, which may have different values 
at different occasions. 

We begin with the following well known fact [2] that the Cahn-Hilliard equation 
) can be interpreted as the ^-gradient flow for the Cahn-Hilliard energy 
functional 

( 2 . 1 ) Mv):= + 


(1.1)-(1.5 


The following assumptions on the initial datum uq were made in m, they were 
used to derive a priori estimates for the solution of problem ( 1 . 11 -( 1.5). 

General Assumption (GA) 

( 1 ) Assume that mg G (—1,1) where 


( 2 . 2 ) 


mo := / uo{x)dx. 

I^‘l Ja 


(2) There exists a nonnegative constant ui such that 


(2.3) 




(3) There exists nonnegative constants cr 2 j erg and 0-4 such that 


(2.4) 


-eAuo-be /(Mo)|| 4 ^qQ) < C'e -^ = 0,1,2. 


Under the above assumptions, the following solution estimates were proved in 


Proposition 2.1. The solution u of problem (1.1)-(1.5) satisfies the following 
energy estimates: 


(2.5) 

( 2 . 6 ) 

(2.7) 


ess sup (^|lVu||i 2 -b -||P(u)||li) 

tG[o,oo) e ^ 

ess sup ||m||^4 < C{1 + Jeiuo)), 

[0,cxd) 

ess sup \\u^ — 1||^2 < CeJe{uo). 

tG [0,cxd) 


ir \\Ms)\\h-i ds 
ir Il^“'(■s)lli 2 ds 


— *^e(lIo); 


Moreover, suppose that {2 


then u satisfies 
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the additional estimates: 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

(2.14) 

(2.15) 

(2.16) 


/ u{x,t)dx = mo Vi > 0 , 

Jn 

POO 

^0 

pOO 

^0 


ess sup \\utr^_^ 

*e[0.oo) ^ / ||^^^||2^^^<^^-max{2.i+3.2.3}^ 

ess sup ||Vy;||^2 Jo 

ie[0,oo) 


ess sup ||AM||i2 < + 

ie[o,cxD) 

ess sup ||VAw||i2 < 

ie[o,cxD) 


!! A ^ ^^-max{2., + |,2.3+i.2.2 + l}^ 

iJo tG[0,oo) 

pOO 

e / WAutWhds + ess sup WutWh < Ce-”""l2ai + f ,2.3+i.2<T3+4.2a4}^ 
^0 te[o,oo) 

pOO 

/ ||A-^Mtt(s)||^-i(is < C'e-“ax{10ai+10,4<Ti+2a3+5.2c73-l}^ 

^0 


Furthermore, if there exists 0-5 > 0 such that 


(2.17) 


lim II VMt(s) 11^2 <Ce 

S-J. 0 + 


then there hold for d = 2, 3, 


(2.18) 

(2.19) 

( 2 . 20 ) 


ess sup II Vitt||^2 + e 

ie[0,oo) 



||VAMt||^2ds < Cpo{e,d), 



\\utt\\H-ids < Cpi{e,d), 


ess sup ||A^u||i2 < Cp2{e), 

ie[0,oo) 


where 

._ g-g^ max{2cri+5,2o.3+2}-max{2CTi + ^,2(T3+|,20-2+4} _j_ g-2o5 
_|_ g— max{2oi+7,2o3+4} 

pi(e,d) := epo{e,d), 

n (e) ■= £“ '^3+1,o-2+f,04 + 1} 


The next lemma concerns with a lower bound estimate for the principal eigenvalue 
of the linearized Cahn-Hilliard operator, a proof of this lemma can be found in [ 5 ]. 

Lemma 2.2. Suppose that (2.2)-(2.4) hold. Given a smooth initial curve/surface 
To, let uo he a smooth function satisfying Tq = {a; S fl;uo{x) = 0} and some profile 


5 




described in m- Let u be the solution to problem jnii-drst- Define Cch as 


( 2 . 21 ) 


Cch :=A(eA-^/'(u)/). 


Then there exists 0 < eo << 1 cind a positive constant Cq such that the principle 
eigenvalue of the linearized Cahn-Hilliard operator Cch satisfies 


( 2 . 22 ) 

for t S [0, T] and e S (0, eo) • 


^ e||VV'lli2 + ^ ^ 

XcH ■= inf —— 1,^ 11^ —- > -Co 

o^tpeH^iii) ^ 

Aw—ijj 




L2 


Remark 1. (a) A discrete generalization of (2.22) on (7° finite element spaces 
was proved in Mm- It plays a pivotal role in the nonstandard convergence analysis 
of m m- In the next section, we shall prove another discrete generalization of 
(2.12) on the DG finite element space. 


(b) The restriction on the initial function uq is needed to guarantee that the so¬ 
lution u{t) satisfies certain profile at later time t > 0 which is required in the proof 
o/ W- One example of admissible initial functions is uq = tanh C^"*-^^ ), where do{x) 
stands for the signed distance function to the initial interface Fq. Such a uq is smooth 
when Fq is smooth. 

Next lemma can be regarded as a nonlinear generalization of the classical discrete 
Gronwall lemma. It gives an upper bound estimate for a discrete sequence which 
satisfies a nonlinear inequality with Bernoulli-type nonlinearity, which will be utilized 
crucially in the next section. A proof of this lemma can be found in [18] and its 
differential counterpart can be seen in m- 

Lemma 2.3. Let {Si\c>i he a positive nondecreasing sequence, {be}£>i and 
{ki}i>i be nonnegative sequences, and p > 1 be a constant. If 


(2.23) 

(2.24) 
then 

(2.25) 
where 

(2.26) 


Si+i - Si < biSi + keSf for I 


i-\ 


P-b (1 - p) ^ A:sa^+i > 0 fori >2, 




*5^ < — <! *51 P + {l-p)Y^ fesOs+i 

I S=1 




for i >2, 


Oi := 


e-i 

nr^- 


3. Fully discrete MIP-DG approximations. In this section we present and 


analyze two fully discrete MIP-DG methods for the Cahn-Hilliard problem (1.1)-(1.5). 


The primary goal of this section is to derive error estimates for the DG solutions that 
depend on only in low polynomial orders, instead of exponential orders. As in the 
finite element case (cf. [H]), the crux is to establish a discrete spectrum estimate for 
the linearized Cahn-Hilliard operator on the DG space. 
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3.1. Formulations of the MIP-DG method. Let Th = {K}k£Q, be a quasi- 
uniform triangulation of 12 parameterized by /i > 0. For any triangle/tetrahedron 
K € Th, we define to be the diameter of iL, and h := max^gj-^ hx- The standard 
broken Sobolev space is defined as 

(3.1) H^Th) ■■= {u e l2(12); VK e %, v\k S H^K)}. 

For any K G Th, Pr{K) denotes the set of all polynomials of degree at most r(> 1) 
on the element K, and the DG finite element space 14, is defined as 


(3.2) Vh := {v e L"(r2); VX S Th, v\k G Pr{K)). 

Let Lq denote the set of functions in L^(r2) with zero mean, and let 14 := 14nLg. 
We also define to be the set of all interior edges/faces oi Th, Eh t)e the set of 
all boundary edges/faces of 7/i on F = dEl, and Eh ■= E^ U Ej^. Let e be an interior 
edge shared by two elements Ki and K 2 . For a scalar function v, define 

{v} = ^{v\k + v\k'), [v]=v\k-v\k', oneGE^, 


where K is Ki or K 2 , whichever has the bigger global labeling and K' is the other. 
The L^-inner product for piecewise functions over the mesh Th is naturally defined by 


{'v,w)rH 



vwdx. 


Let 0 < 4 < 2^1 < • • • < tju = T be a partition of the interval [0, T] with time 
step k = tn+i — tn- Our fully discrete MIP-DG methods are defined as follows: for 
any 1 < m < M, (G™, IF™) G Vh x Vh are given by 

(3.3) (dtC/™, 77 ) + a„(lF™,r/) =0 V 77 G F?., 

(3.4) ,ah{U”^,v)+^{r,v)-iW^,v)=0 WvGVn, 

e 

where 

(3.5) ah{u,v) = / Tu ■ Tv dx — {Tu-ne}[v]ds 

KeTh 

f{Tvne}[u]ds+ ( '^[u][v]ds, 
e&ei eeei ^ 

and CTg > 0 is the penalty parameter. There are two choices of /™ considered in this 
paper, namely 

/™ = (C7™)3 - C7™-i and /™ = (G™)^ - C/™, 


which lead to the energy-splitting scheme and fully implicit scheme respectively, dt 
is the (backward) difference operator defined by dfU"^ := (C/™ — U'^~^)/k and := 
PhUQ (or QhUo) is the starting value, with the finite element (or L^) projection 
Ph (or Qh) to be defined below. We refer to [13] for a discussion why a continuous 
projection is needed for the initial condition. We remark that only the fully implicit 
case was considered in [UlTi] for the mixed finite element method. 
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In order to analyze the stability of (3.3)-(3.4|, we need some preparations. First, 
we introduce three proiection operators that will be needed to derive the error es- 

Ph : 


timates in section 3.4 
defined by 


Vh denotes the elliptic projection operator 


(3.6) ah{u-PhU,Vh) + {u-PhU,Vh) =0 G 14, 

which has the following approximation properties (see 0): 

(3.7) \\V - PHvhnn) + h\\w{v - Phv)\\mn) < 

Here r := min{l,r} — min{l, r — 1}. 

Let Ph : H^{Th) —t Sh ■= 14 H (^^(H) denote the standard continuous finite 
element elliptic projection, which is the counterpart of projection Ph- It has the 
following well-known property [TilfT^] : 

(3.9) \\u-Phu\\L^ <Ch^-^\u\\H^. 


Next, for any DG function G 14, we define its continuous finite element 
projection 4'^'® G by 

Uhi^h^.Vh) = ahi^h,Vh) yvh€Sh, 


(3.10) 

where 


dhiu, v) = tthiu, v) -I- a(u, v), 


and a is a parameter that will be specified later in section 3.3 

A mesh-dependent H~^ norm will also be needed. To the end, we introduce 
the inverse discrete Laplace operator ■ Vh ^ Vh as follows: given ( G 14, let 
Aj^^^ G Vh such that 

(3.11) ah{-A~^(,Wh) = (C,Wh) ywhGVh- 

We note that is well defined provided that (7° > cr° for some positive number 
(T° and for all e G £h because this condition ensures the coercivity of the DG bilinear 
form a/i(-, •). 

We then define “-1” inner product by 

(3.12) ■■= ah{-A-\,-A-^0 = 

and the induced mesh-dependent H~^ norm is given by 

(C,0 


(3.13) 

where 


IlCII-i.h := \j{C,0-i,h = sup 


•= \/aft(4 0- The following properties can be easily verified (cf. [3]): 

(3.14) |(C,C)l<IICII-i,/.llieilla VCG14, C€l4, 

(3.15) ||CII-1./.<C'||CIU^ VCG14, 

and, if Th is quasi-uniform, then 

(3.16) IICIU 2 


VCg 14. 










3.2. Discrete energy law and well-posedness. In this subsection we first 
establish a discrete energy law, which mimics the differential energy law, for both 


fully discrete MIP-DG methods defined in (3.3)-(3.4). Based on this discrete energy 


law, we prove the existence and uniqueness of solutions to the MIP-DG methods by 
recasting the schemes as convex minimization problems at each time step. It turns 
out that the energy-splitting scheme is unconditionally stable but the fully implicit 
scheme is only conditionally stable. 

Theorem 3.1. Let G Vh x Vh be a solution to scheme (3.3)-(3.4). 

The following energy law holds for any h,k > 0 : 

(3.17) + k WdtUnWh + fc" E I l\\\dtU"^\\\l + 

tn—l m—1 K. 

+ ± = E,{u°) 


for all 1 < £ < M, where 


(3.18) 


En{U) :=l||[/2-l||i. + | 


f 


Note that the sign “±” in ([3T7| takes “-k” when = {U'^f and when 
= - U^. 


The proof of the above theorem follows from taking rj = —A~^dtU"^ in (3.3 1 and 
V = dtU™' in ( |3.4[ ), adding the resulting two equations and combining like terms. We 
leave the detailed calculations to the interested reader. 

Corollary 3.2. Let > 0 be a sufficiently large constant. Suppose that 
o'e > crj for all e G Sfi- Then scheme (3.3)-(3.41 is stable for all h,k > 0 when 
~ U'^~^ and is stable for h > 0 and k = 0(e®) when = {U^Y ~ ■ 


Proof. The first case holds trivially from (3.17). In the second case, the “bad 
term” \\dtU"^\\L 2 can be controlled by the “good terms” ^ and IIIG^IIIa by 


using the norm interpolation inequality (3.40) provided that k = O(e^). □ 

Theorem 3.3. Suppose that for all e G £h- Then scheme (3.3)-(3.4| 

has a unique solution (t/™,IT'") at each time step for for all h,k > 0 in the case 
/”* = {U^Y ~ and for h > 0 and k = O(e^) in the case Z™ = {IT^Y ~ . 


Proof. Setting p = — in 


we get 


{d,U-,v)_^.+{W-,v)=0. 


Adding the above equation to (|3.4[) yields 


{dtU^,v)_^^ + ean{U^,v) + -(/™,u) = 0. 


Hence, G’" satisfies 


(3.19) 


(C/™,u)_^, + fcea4t/-,u) + -ir^v) = {U-\v) 


In the case /"* = {U'^Y ~ it is easy to check that (3.19) can be recast 

as a convex minimization problem (cf. [31 I13j i whose well-posedness holds for all 
h,k > 0. Hence, in this case there is a unique solution 17™ to (3.3)-(3.4). On the 
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other hand, when /"* = — t/™, there is an extra term —ke comes 

out from the nonlinear term in (3.19). This extra term contributes a “bad term” 
to the functional of the minimization problem. Again, this term can 
be controlled by the “good terms” and |||17'"|||a in the functional by using 

the norm interpolation inequality (3.40), provided that k = O(e^). Hence, in the case 
/”* = {U^Y — U™, there is a unique solution 17™ to (3.3)-(3.4) for all > 0 and 
k = O(e^). The proof is complete. □ 

3.3. Discrete spectrum estimate on the DG space. In this subsection, we 
shall establish a discrete spectrum estimate for the linearized Cahn-Hilliard operator 
on the DG space, which plays a vital role in our error estimates. 

To the end, we first state a slightly modified version of a discrete spectrum esti¬ 
mate for the linearized Cahn-Hilliard operator on the continuous finite element space 
first proved in [Hin]. Due to the close similarity, we omit the proof of this modified 
version and refer the interested reader to [II US]. 

Lemma 3.4. Suppose the assumptions of Lemma 2.2 hold, and Cq is the same as 


in ( 2 . 22 ). Cl and C 2 are defined by 




(3.20) 

(3.21) ||u - P.u|U=o((o,t);l~) < 

Then there exists 0 < ei << 1 such that, for any e G (0, ei), there holds 


(3.22) 


Ag| = inf 


\yA-^Yh\\h 


^ “(C'o + 1 ); 


provided that h satisfies 
(3.23) 

Here : Lq{LI) —> H^(J1) H To(^) denotes the inverse Laplace operator. 

We are now ready to state the discrete spectrum estimate on the DG space. 
Proposition 3.5. Suppose the assumptions of Lemma \2.^ hold. Let u he the 
solution of (!.!)-( 1.5) and P^u denote its DG elliptic projection. Assume 


(3.24) 


ess sup ||u||wi+^.“ < C'e 

i:G[0,CXD) 


for a constant 7 , then there exists 0 < £2 << 1 ctnd an e-independent and h- 
independent constant cq > 0, such that for any e G (0, £2), there holds 


Agg = inf 


eohi^h, ^h) + ^-^ifiPhu)^h, ^h) 
l|VA-i<l>,|| 2 , 


(3.25) 

05^<i.hGi2(n)nVh 
provided that h satisfies the constraints 

(3.26) < (C'iC'2)”^£“’"’'^‘"' + ^’‘"=^+^g 

(3.27) In hf < {CiC3)-Y^+^. 

where Ci and C 2 are same as in Lemma \T^ f and C 3 are defined by 
f = min{l, r} — min{l, r — 1}, 

||u - P^.ir||L~((o.T);L~) < In hfe-P 


> -Co, 
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Proof. By Proposition 2 in |T4j, under the mesh constraint (3.26), we have 

(3.28) WfiPhu) - < e3. 


Similarly, under the mesh condition (3.27), we can show that for any e > 0, there 
holds 

(3.29) Wf'iPhu) - /'(w)|U=.((o, t);L») < 

It follows from ( 3.28| ) and (3.29) that 

(3.30) |1/'(Am) -/'(Au)||l~((o.t);L-) < 2e^ and f {Phu) > f {Phu) - 2f?. 
Therefore, 

1 _ g3 

(3.31) tah{<^h, $/.) +- {f{Phu)<^h, <^h) 


P ‘h/i) + 

l-e3 


1-£3 

e 


(/'(P,u)$^,$^)-2e2(l-e3)||ci,^ll^^ 


= e- 


1 _ £1 
^ 2 


a/*($h,$/j) + — -{f{Phu)^h,^h) 


- 26^(1-e3)||$.||i.+ 
l-e3 

-ahi^h, ‘h/t) + 


1 - 2 
l-e3 


2-e3 

l-e3 

\FE\ 2 , 


dhi^h, ‘h/t) 

/ fiPhu)(^{<^>hf - dx 


nPnu)i<^^^rdx - 2e^il - e3)|i$,||2, + $,). 

I Z — 6 


Next, we derive a lower bound for each of the first two terms on the right-hand 
side of (3.31). Notice that the first term can be rewritten as 

(3.32) $^) = - $^^) + 2a^($,, $^^) - 

= a.(d>. - <i>r, ‘i’/. - + l|V$nii. + 2a\\<^r - <i>h\\h 

+ 2a{<^>f^-^h,^h)- 

To bound from above, we consider the following auxiliary problem; 

M^,x) = i^h-^r,x) Vxepi(ii). 

For cTg > a° for all e G £h, the above problem has a unique solution cj) e for 

0 < 0 < 1 such that 

(3.33) \mm+o(n)<C\\^h-^i^\\L^ for 0G (0,1]. 

By the definition of we immediately get the following Galerkin orthogonality: 

dh{<^h-,Xh) yXh&Sh, 

It follows from the duality argument (cf. [2T1 Theorem 2.14]) that 

(3.34) 11$;, - $rili2 < 

< - <hr, 
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For all h satisfying Clr^a < 1, we get 

(^h'^0 

(3.35) 11$, - - <^). 


Now the last term on the right-hand side of (3.32) can be bounded as follows: 

(3.36) 2a($f^> -2a\\^l^ - ^hWl-W^hWl- 


> -2a\ 


1 - CK^^a 

hFE ^FE\ 




> -^ahi^h - $/» - $C^) - i^c!^2ej ‘^h\\l2. 


The second term on the right-hand side of (3.31) can be bounded by 

(3.37) [ /'(P,w)(($,)2-($H')rfa:>-C /|($,)2-($n"|da: 

Jq J q, 

= -C ( -($,-$^^)%2$,($,-$^^) dx 
Jn 


>-cii$,-$rili2- 


3(1 - £3) 

1_ £ 

2 


\Mh-c- 


1 - 


£ 3(1 - £ 3 ) 

Here we have used the facts that 

(3.38) ||u|U=o((o.t);l~) < C, \f{Phu)\ < \fiu)\ + £3 < C. 

Substituting ( 3.35| ) into ( |3.37| ) yields 
1-£3 




(3.39) [ /'(P,u)(($,)2 - ($^^)2) dx 

e Jo. 


> -73 


.(1-£3) 

1_ Jd 
^ 2 




^(1-£3) 

1 - 
^ 2 


W^hWl., 


where 


£3a-"£3))’ 


£ 3(1 -£ 3 ), 


. ^ gr ^ ^ 

^ 1 — Ch^^a £(1 — £ 3 ) 

and h is chosen small enough such that 73 < 1/4. 

The term ||$/t |||2 can be bounded by 

(3.40) |l$,||i 2 = ($,,$,) = a,(A/3$^^$^) < a,(A/^$,, A/^$,)5a,($,,$,)5 

^ + ^ah{^h,^h) 

for any constant p > 0 . 

Adding the fifth term on the right-hand side of (3.31), the last term on the 
right-hand side of (3.36) and that of (3.39), we get for all h satisfying 2Ca'^h‘^^/{I — 
CK^^a) < £ 


(3.41) - 


;(1 - £ 3 ) 2Ca^h^^ 3e^{l - e^) 


1 _ <1 1 - Ch^<^a 


1 - 




1 _ <7 
2 


> -- 


2 ( 2 -£3) 


^ J - Cah{A^ 1^ 
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Combining (3.32), (3.36), (3.39) and (3.41) with (3.31), we have 


(3.42) eahi<^h,^h) + 


1 -e^ 


f {Phu){^hf dx 


> 


+ 


e(l — e^) __ ^ ^FE\ I 2ae(l — e^) 


4- 2e3 

^3' 




:3 
(3 

\\V7,f.FEn2 


1 - 




1 -' 
l-e3 




r{p,u){<i>i^fdx+ 


2 (2-e3) 




Applying the spectrum estimate (3.22), we get 
l-e3 


1 - 


iiv$nii. + — [ nPHu){<i>rf dx 


^ ' iiv$nii. + [ npHux^^^fdx 


i-f 


In 


^ 2 


which together with (3.42) implies that 

l_e3 , 


(3.43) eahi^h,^h) + 


f'{Phu){^h)dx 


> -CaniA-^^^, A-ich;,) - C|| VA-^$rilL^ + 


-i^FE \\2 I 2ae(l e^) 


1 _ SI 
2 




By the stability of A we have 

iivA-i(ci>, - <^)iii, < aii$, - 

which together with the triangle inequality yields 

||VA-i<l>nii. < 2||VA-i<l>;,||i. + 2^11$, - 

Similarly, since A^^^h is the elliptic projection of A~^$?i, there holds 


Therefore, choosing a = 0{Ce ^), (3.43) can be further reduced into 

eah($/t, $/i) + ^^ [ fiPhu)i<^hP dx > -Co|| VA "^$,,|||2 

e Jn 


for some cq > 0. This proves (3.25), and the proof is complete. □ 
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3.4. Error analysis. In this subsection, we shall derive some optimal error 


estimates for the proposed MIP-DG schemes (3.3|-(3.4|, in which the constants in 
the error bounds depend on e~^ only in low polynomial orders, instead of exponential 
orders. The key to obtaining such refined error bounds is to use the discrete spectrum 


estimate (3.251. In addition, the nonlinear Gronwall inequality presented in Lemma 


|2.3| also plays an important role in the proof. To ease the presentation, we set r = 1 
in this subsection and section and generalization to r > 1 can be proven similarly. 

The main results of this subsection are stated in the following theorem. 
Theorem 3.6. Let {(G™, IT™)}^^q be the solution of scheme {3.S)-{3.4) with 
r = 1. Suppose that (GA) holds and > cr° for all e G Sh, and define 


(3.44) 

(3.45) 


■= g-max{2CTi + ^,2(T3+| ,2(72+4,2(T4}-4 

r{h, k] e, d, a^) := /c^pi(e; d) + h^p 3 {e). 


Then, under the following mesh and starting value conditions: 


(3.46) 

(3.47) 

(3.48) 

(3.49) 

(3.50) 

(3.51) 

(3.52) 


h^+nin hf < {CiC3)-h^+^, 

k<€^ when /“ = {U'^ f - 

+»< 

- 8-4e3 ’ 

(t7°,l) = (uo,l), 

Iko — < C'^^||mo||h 2 , 


there hold the error estimates 


M 


(3.53) max Mu) - U"^\\-i,h + ( k’^WdtMU) - UM\-i,h) ' 

0<m<M \ / 

m—\ 

< Cvih^ k; e, d, cr^) 2 , 

M 1 

(3.54) {kJ2Mtm)-U^My 

m—1 

< + e-‘^r{h,k;e,d,a,)y, 

M 1 

(3.55) {kY,M{u{U)-U^)My 

m—1 

Moreover, if the starting value M satisfies 

||mo — U°\\l 2 < Ch‘^\\uo\\H^, 


(3.56) 
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then there hold 


M 1 

(3.57) max \\u{tjn) - U^\\l 2 + (kk\\dt{u{tm) - U^)\\l2)'‘ 

0<m<M \ / 

ra—\ 

k ^ 1 

771=1 

< C(h'^p3{e)^ + e~ir{h,k;€,d,a,)i^ 

(3.58) \\u{tm) - f7'"||i=o 

0<m<M 

< In/ije”''' + h~^e~ir{h, k;€,d,ai)^^ . 

Furthermore, suppose that the starting value satisfies 

(3.59) \\PhWo-W°\\L2<Ch^ 


for some f > 1, and there exists a constant 7 ' such that 


(3.60) 


ess sup ||w||^y 2 ,oo <Ce'^ 

iG[0,C5o) 


then we have 

(3.61) max \\w(tm) - < c(h'^p 3 {e) + h^ 

0<m<M \ 

+ k~^e~^r{h, k;e,d,ai)^^, 

(3.62) max \\w{tm) — < c(h~^ (k~^e~^r{h,k;e,d,ai)^ + h^\ 

0<m<M V V / 

+ h'^\\nh\e-^'y 


Proof. In the following, we only give a proof for the convex splitting scheme 
corresponding to /"* = (tt™)^ — in (3.13) because the proof for the fully implicit 
scheme with /”* = (u™)^ — is almost same. Since the proof is long, we divide it 
into four steps. 

Step 1: It is obvious that equations (fTIt-dOt imply that 

(3.63) {utftm), Vh) + ah{w{tm), Vh) =0 & 14, 


(3.64) 


eah{u{t 

m ),Vh) + -{f{u{tm)),Vh) = {w{tm),Vh) VUft G 14. 


Define error functions E"^ := u{tm) — and G"* := w{tm) — ID™. Subtracting (3.3) 
from (3.63) and (3.4) from ( |3.64[ ) yield the following error equations: 

(3.65) {dtE^,ph) + ah{G'^,ph) = {R{utt,rn),ph) V???, G 14, 

1 


(3.66) ea^(i?™, vu) + - {f{u{t^)) - f(U^), vu) = (G™, vu) 

where 


V'C/i G 14, 


Rfuturn) := 


(s - tm-i)utt[s) ds. 
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It follows from (2.191 that 


k'^\\R{utum)\\jj-i < -'^U {s - dsj U \\utt{s)\\H-i dsj 

< C'fcVi(e, d)- 

Introduce the error decompositions 
(3.67) = 0 ™ + C?”" = A"" + 

where 


0 ™ := := 

:= d/™ := - Vh™. 

Using the definition of the operator Ph in ( |3.6[ ), ( |3.65[ )-( [3.66[ ) can be rewritten as 

(3.68) {dt<^"^,r]i,)+ah{'i’"",Vh) =-{dt'd"",Vh) + {R{utt,m),r]h) V% G U/,, 

(3.69) eaft,($’", vt) + - {f{u{tm)) - /™, vt) = (5''", Vh) + (A™, vh) ^vh G 14. 


Setting r]h = — in (3.68) and Vh = d)™ in (3.69), adding the resulting 

equations and summing over m from 1 to £, we get 

e 

(3.70) A^i$^) + ^ - A;)!^'"-!) 


+ 2kJ2 ch™) + 2kJ2- ifHtm)) - r, < 5 ”^) 


m—1 


m—1 


= 2fc ^ [{R{uu, m), - {dt&^,+ (A™, d>™)) 

m—1 

+ a,(A-i$°,A-ici>°). 


Step 2: For ( 7 ° > tr® for all e G £h, the first long term on the right-hand side of 


(3.70) can be bounded as follows 


(3.71) 2k Y, {{R{uu, m), + {dtQ^, + (A™, d>™)) 

m—1 

i 

<CkY (||P(««;m)||^_. + ||d*0-||^_. + (1 - e3)e-4||A™||^_.) 

m—1 

t 4 

+ kY $-)) 

m=l 

i 4 

< fc ^ (a^(A-i<i>™, A-i$™) + ^^a„(<l>™, $“)) 

m—1 

+ C'^fc^pi(e, d) -b p3{e)^ , 
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where we have used (2.15) and the following facts [T(7| : 


||m - Phu\\H~i < Ch^\\u\\H^, llw - Phw\\H-i < Ch^\\w\\H^. 


We now bound the last term on the left-hand side of (3.70). By the definition of 
/™, we have 


- /"* = f{u{tm)) - f{Phu{tm)) + f {Phu{t^)) - 

> -\f{u{t^)) - f{Phu{tm)) I + {Phu{tm)f - Phu{tm) " (C/™)" + 

> -Cl©'" I + ({PhU{t^)f + PhUit^) - $» _ kdtU^ 

> -C|0'"| + r{Phu{tm)) - 3Phu{trn) - kdtU^. 


By the discrete energy law (3.17), (3.13) and (3.40), we obtain for any 1 < £ < M 


(3.72) 2fc^ i(/(u(U) 

m—1 

I I 

m—1 m—1 

lit 

- ^ E ^ ■ T ^ iid>™iu 

m—1 m—1 m—1 

I It 

> 2fc E \{j'{Py.u(tm))A‘^n^) + ^ E - V ^ 

771=1 771=1 771 = 1 

- E 

m—1 

t it 

> 2fc E \(!\PKu{t^))A^n'^) + Y E - V ^ 

771=1 771=1 771 = 1 

4 t 








Substituting (3.711 and (3.72) into (3.70) we get 


(3.73) ^ 




l-e^ 
e 




m—l 


m—1 


<CkY, ahiA-^^-^, A-^‘^>^) + — E 

m—l m—l 

£ 

- lOfce" ^ {f{Phu{t^))<^”^, $-) + C(fc"pi(e; d) + h^p^{e)) 

m—l 

+ c(ji^e ®||M||L2((o,T);i/=(n) + ^Eh{U°)j. 


Step 3: To control the second term on the right-hand side of (3.73), we appeal to 
the following Gagliardo-Nirenberg inequality [1]: 


^Ili3(ic) < C(|1V^;|| 


Thus we get 
(3.74) 


LHK)r\\L^[K) + \\^ 


I 


Ili2(if)) 


V7f G Th. 


^ E E w^^'^winn) + V E 

m—l m—l m—l 


4(l + rf) V^ll^' 

-I- Ce i-'i k ^ ||4> 


2(6-d) 
m II 4—d 
L 2 


< 


Cfc 


A-e3 




4( 1 + d) ^> 11 

+ Ce-^^k^\\<^' 


2(6-d) 
m II 4-d 
L 2 


m—l 


(3.75) 


The third item on the right-hand side of (3.73) can be bounded by 


- + kCahiA-^<S>"^, A-'$™). 


Again, here we have used (3.40). 

Finally, for the third term on the left-hand side of (3.73), we utilize the discrete 
spectrum estimate (3.25) to bound it from below as follows: 

(3.76) ea44>™,4>™) + ^^(/'(P,M(t^))$’",$’") > -co||VA-i4>™||2,. 
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By the stability of A ^ and (3.40), we also have 


(3.77) co||VA-i$-||2, <c|l$-||2, < 


Step 4'- Substituting (3.74), (3.75), (3.76), (3.77) into (3.73), we get 


(3.78) 


f. 

a;,(A-i4>^ A-i$0 + ^ a^(A-i$'" - A-i4>”^-\ A-i$- _ 

m—l 

+ ^ ^ ^ ||ci>-||4, 

ra—1 m—l 


f / 

Ck 

<CkY^ a,(A-i$-, A-i$-) + — E 

m—l m—l 

e 

4(l + d') 2(6-d) 

+ Ce- —E d) + hV3(e)) 

771=1 

+ C^h®e ®||M||| 2 ((o,T);ff=(n) + ^Eh{U°)^. 


By discrete energy law (3.17), General Assumption (2.3), stability of ellip¬ 
tic projection, L°° stability(or L°° error estimate and triangle inequality) of elliptic 
projection, we can get for any 0 < i < M 


£ 

\\U‘\\l2 <kJ2 \\dtUn\L- + \\U°\\l^ < 

771=1 

Since the projection of u is bounded, then for any 0 < £ < M 


(3.79) |1$^||l 2 < C'e"‘"b 

We point out that the exponent for ||<1>'"|L2 is which is bigger than 3 for 

d = 2, 3. By we have 




I"* < 

1 ^2—'-^'^ II ^ II L2 ’ 




lL 2 - 


< I Id)™ I 


Lo¬ 


using the Schwarz and Young’s inequalities, we have 

3 

(3.80) =a^(-A;^i$’",$’")i 

<a;,(A-i$”^,A-i$’")ta„($™,$™)i 

< gT^+<^i+2(d-2)cri_e_a/i(<i)'", $”") 

1 — 
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Therefore, (3.78) becomes 


(3.81) ah{A-^<^>\ - A-i$™-i) 


m—1 


m—l m—1 

i 

<C'fc^a/,(A-i<i>™,A-i<i>™) 

771=1 

m—1 

+ C'(fc^pi(e; d) +/i®p3(e)) + C^/i®e ^||M||^2((o,T);ff»(n) + ^Eh{U^ 
i 

<Cfc^a;,(A-i<i>™,A-i<i>™) 


+ C'fce-^^-2<Ti-4(d-2),Ti J2 a,,(A^i$’",A;^i$’")^ 


771=1 


+ C{k'^pi{e; d) + h^psie)). 

On noting that C/™ can be written as 

i 

(3.82) = kY^ dtU”^ + U°, 

m—1 


then by (2.3) and (3.17), we get 


(3.83) \\U^\\-iM <kJ2 \\dtUn\-i.h + ||C/°||-i,„ < Ce' 

771=1 

Using the boundedness of the projection, we have 

(3.84) 


Also, (3.81) can be written in the following equivalent form 


(3.85) ahiA-;;^^^ A^i$^) + ^ ahiA^^<^^ - 

771 = 1 

4^ ^ OjL ^ 

+ E $-) + - ^ ||$-||i 4 < Ml + M 2 , 

771=1 771=1 


20 










where 


z-\ 

(3.86) Ml := Ck ^ 

m—\ 

i-1 

m—1 

+ C{k'^pi{e-, d) + h^p 3 {e)), 

(3.87) M 2 := Cfca„(A)(l$^ A))i$^) 

It is easy to check that 

(3.88) M 2 < provided that fc < 

Under this restriction, we have 


i 

(3.89) a,,(A))i$^ A^i$^) + 2 ^ aft(A;(i$™ - A;(i$™ - A)(i$™-i) 

m—X 

+ E <i>™) + ^ ^ ||<i>™||i4 

m=l m=l 

£-1 

< 2Ck ^ oi,(A^i$'", A)^!^’”) + 2C{k^p,ie; d) + h^^p^ie)) 

m—1 

m=l 

^-1 

<CkY, A)(i<I>™) + C(fc2pi(e; d) + h^p^{e)) 

m—1 

m=l 

Define the slack variable dg > 0 such that 


e 

(3.90) ahiA-^^\ A-i$^) + 2 ^ ah{A-^^^ - A-^<^^-\ A-^^^ - A-^^^-i) 

m—1 

+ ^ E $-) + ^ ^ ||<I>-||i4) + d, 

m=l m—1 

i-1 

= CkY A)(i<I>™) + C'(fc2pi(e; d) + /iV3(e)) 

m—1 

i-1 

+ C'fce-^)S^-2'-i-4(d-2)<Ti ^ ai,(A“i$'",A-i$’”)3. 

m—1 
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We also define by 


(3.91) St = di + 2Y^ 

m—1 

+ a^(A;^i$^ a;) 1$«) + ^ a,(<i>™,Ci>™) + ^ ^ ||$” 

m=l m=l 


lli^, 


and equation (3.90) shows that 

Si=C{k‘^pi{e-d) + h^p^{e)). 

Then 

(3.92) Si+x-+ > 1. 


Applying Lemma 2.3 to defined above, we obtain \/^ > 1, 

(3.93) St < S'f^ - 


provided that 

(3.94) 




£-1 


We note that Os (1 < s < £) are all bounded as fc —)■ 0, therefore, (3.94) holds under 
the mesh constraint stated in the theorem. It follows from ( |3.51 ) and (3.52) that 

(3.95) St < 2aJ^Si < C{k^pi{e;d) + h^psie)). 


Then (3.53) follows from the triangle inequality on if"* = 0”* + $"*. (3.55) is 
obtained by taking the test function ph = 4>"* in (3.68) and Vh = 4>”* in (3.69), and 
(3.54) is a consequence of the Poincare inequality. 

Now setting ph = ‘h™ in (3.68) and Vh = —^4'”* in (3.69), and adding the resulting 
equations yield 

(3.96) IdtW^nih + ^Wdt^nih + h'i’nih = - /(t/’"),4/"*) 

+ {R{uu; m), $"*) - (dt©”*, $”*) - i (A”*, 4-"*). 

The last three terms on the right-hand side of ( 3.96[ ) can be bounded in the same way 
as in (3.71), and the first term can be controlled as 


(3.97) 




<-\\'i>^\\h + ^\\E^\\h. 


1 

Ye' 


Multiplying both sides of (3.96) by k and summing over m from 1 to M yield the 
desired estimate (3.57). Estimate (3.58) follows from an applications of the following 
inverse inequality: 

(3.98) 11$" 
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and the following L°° estimate for the elliptic projection: 
(3.99) \\u - PhuWl-- < Ch^\\nh\\\u\\w-^-- Vu S 


Finally, it is well known that there holds the following estimate for the elliptic 
projection operator: 


(3.100) 


0<m<M 

Using the identity 
(3.101) 
we get 

1 


M 


max IjA^IU. + fc E k\\dtA-^\\h < Ch^P 2 {e). 


m—0 




(3.102) 




M ||2 

Il2 


M , M .. 


\l^ 


m—l 

M 


< A: E ( jimnih + 


m—l 


The first term on the right hand side of (3.102) can be absorbed by the second term 


on the left hand side of (3.102). The second tern on the right hand side of (3.102) 


has been obtained in (3.57). Estimate (3.61) for VF"* then follows from (3.100) and 


(3.102). (3.62) follows from an application of the triangle inequality, the inverse 


inequality, and (3.99). This completes the proof. □ 


4. Convergence of numerical interfaces. In this section, we prove that the 
numerical interface defined as the zero level set of the hnite element interpolation of 
the solution converges to the moving interface of the Hele-Shaw problem under 
the assumption that the Hele-Shaw problem has a unique global (in time) classical 
solution. To the end, we first cite the following PDE convergence result proved in [5]. 

Theorem 4.1. Let LI be a given smooth domain and Too be a smooth closed 
hypersurface in LI. Suppose that the Hele-Shaw problem starting from Too has a unique 
smooth solution {w,r := Uo<t<T(^t ^ {^})) ^de time interval [0,T] such thatVt C H 
for all t G [0, T], Then there exists a family of smooth functions {uo}o<e<i which are 
uniformly bounded in e G (0,1] and (x, f) G LIt, such that ifu'^ solves the Cahn-Hilliard 
problem ( 1 . 1 ) , then 

if{x, t) GO 
if{x, t) Gl 

I and O stand for the “inside” and “outside” ofT; 


(i) limu'^(a:,t) = 
£->■0 


uniformly on compact subsets, where 


(ii) hm(e 
£->■0 


-1 


/(m*^) — eAu'^^{x,t) = —w{x,t) uniformly on LIt- 


We note that since t/™ is multi-valued on the edges of the mesh 7)i, its zero- 
level set is not well defined. To avoid this technicality, we use a continuous finite 
element interpolation of t/™ to define the numerical interface. Let U™ G St denote 
the finite element approximation of U™ which is defined using the averaged degrees of 
freedom of U™ as the degrees of freedom for determining 17™ (cf. [II]). The following 
approximation results were proved in Theorem 2.1 of m- 

Theorem 4.2. Let Tk be a conforming mesh consisting of triangles when d = 2, 
and tetrahedra when d = 3. For Vh G Vh, let Vh be the finite element approximation 
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of Vh as defined above. Then for any Vh G Vh and i = 0,1 there holds 


(4.1) 


ik/i—< c'y]] 

KSiTh eG££ 


2 

L 2 (e)> 


where C > 0 is a constant independent of h and Vh but may depend on r and the 
minimal angle Oq of the triangles in Th- 

By the construction, U'^ is expected to be very close to C/™, hence, [/™ should 
also be very close to ujtm )- T his is indeed the case as stated in the following theorem, 
which says that Theorem 


3.6 


also hold for 


Theorem 4.3. LetU'^ denote the solution of scheme (3.1)-(3.14) andU"^ denote 
its finite elem ent approximation as defined above. Then und er the assumptions of 
Theorem 


3.6 the error estimates for 17™ given in Theorem 3.6 are still valid for IT 

( 


in particular, there holds 
(4.2) 


^miix^\\u(tm) - U"^\\L°°{Th.) < Cyh^\\nh\e ^ ^ 


r{h,k-,e,d,ai)^^y 

to see a proof of the 


We omit the proof to save space and refer the reader to 
same nature for the related Allen-Cahn problem. 

We are now ready to state the first main theorem of this section. 

Theorem 4.4. Let {rt}i>o denote the zero level set of the Hele-Shaw problem 
and {U,:^fi,k{x,t),W,:^u,k{xT)) denote the piecewise linear interpolation in time of the 
finite element interpolation {(17™,TT™)} of the DG solution {(U™, W™)}, namely, 


(4.3) 

(4.4) 




t tvn — 1 


U^ix) + 


t-m t ' 




t tm —1 


W™(a:) + 


tm, t 


W'^-^x), 


for tra-i t < tm and 1 < m < M. Then, under the mesh and starting value 
constraints of Theorem 3.6 and k = 0(hf~^) with 7 > 0, we have 

(i) Ue,h,kix,t) 1 uniformly on compact subset ofO, 

(ii) Ue.h,k{x,t) —1 uniformly on compact subset ofT. 

(hi) Moreover, in the case that dimension d = 2, when k = 0{h^), suppose 
that satisfies ||wg — W°||l 2 < Ch^ for some /3 > |, then we have 

We,h,k{x,t) —w{x,t) uniformly on LIt- 
Proof. For any compact set A C O and for any {x, t) G A, we have 


(4.5) 


lU,,^h,k - 1 | < \Ue,h.k - u’'{x,t)\ + |u"(a;,f) - 1 | 

< \Ue,h,k - U%X,t)\L^(^Q.^-j + \u’'{x,t) - 1 |. 


Equation (3.581 of Theorem 
that 


3.6 


infers that there exists a constant 0 < a < such 


(4.6) 


\Ue,h,k - '«"(a;,t)|L=o(nT) < C!h° 


The first term on the right-hand side of (4.51 tends to 0 when e \ 0 (note that 
h,k \ 0, too). The second term converges uniformly to 0 on the compact set A, 
which is ensured by (i) of Theorem 4.1 Hence, the assertion (i) holds. 
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To show (ii), we only need to replace O by X and 1 by —1 in the above proof. 
To prove (ill), under the assumptions k = 0{h^), ( |3.62 ) in Theorem 3.6 implies that 
there exists a positive constant 0 < C < such that 

(4.7) \\W,,h,k-w^\\L^(^nr)<Ch‘^- 
Then by the triangle inequality we obtain for any {x,t) G VIt, 

(4.8) \We,h,kix^i) - (-''")l < |4Te,/i,/c(a:,t) - w"(a:,t)| + |w"(a:,t) - {-w)\, 

< l|W4,/i.fc(a;,i) - ■ic'^(a;,t)||i=o(o^) + \w''{x,t) - {-w)\. 

The first term on the right-hand side of ( |4.8| ) tends to 0 when e \ 0 (note that 
/i, fc 0, too). The second term converges uniformly to 0 in fix, which is ensured by 
(ii) of Theorem 4.1 Thus the assertion (iii) is proved. The proof is complete. □ 

The second main theorem of this section which is given below addresses the con¬ 
vergence of numerical interfaces. 

Theorem 4.5. Let := {x £ LI; Ue^h,k{x,t) = 0} he the zero level set of 

Ue,h,k{x,t), then under the assumptions of Theorem 4-4' we have 


sup dist(a;,rt) 


e \,0 


0 uniformly on [0,T]. 


Proof. For any t] £ (0,1), define the open tabular neighborhood A/), of width 2r] 
of T t as 

(4.9) Mrf := {{x,t) £ LIt; dist(x,rt) < 77 }. 

Let A and B denote the complements of the neighborhood A/), in O and I, respectively, 

i.e. 


A = 0\Mr, and B=I\Mr,. 

Note that A is a compact subset outside Tj and B is a compact subset inside Ft, then 
there exists £3 > 0 , which only depends on 77 , such that for any e G ( 0 , 63 ) 


(4.10) |t/e, 7 i,fe(a;,t) - 1| < 77 V(a;,t)eA, 

(4.11) |t/e,/i,fe(a;,t) + 1| < ?7 y{x,t)&B. 

Now for any t £ [0,T] and x £ from U^^h,k{x,t) = 0 we have 


(4.12) |17,,,,,fc(a;,t)-l| = 1 y{x,t) £ A, 

(4.13) \U,,h,k{x, t) + l\ = l y{x, t) £ B. 


(4.10) and (4.12) imply that (a;,t) is not in A, and ( 4.11| ) and (4.13) imply that {x,f) 
is not in B, then {x,t) must lie in the tubular neighborhood A/);. Therefore, for any 
e G ( 0 , £ 3 ), 


(4.14) 


sup dist(a;,Fj) < 77 uniformly on [0,r]. 


The proof is complete. □ 
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5. Numerical experiments. In this section, we present three two-dimensional 
numerical tests to gauge the performance of the proposed fully discrete MIP-DG 
methods using the linear element (i.e., r = 1). The square domain G = [—1,1]^ is used 
in all three tests and the initial condition is chosen to have the form uq = tanh( ^^^ ), 
where cIq^x) denotes the signed distance from x to the initial interface Tq. 

Our first test uses a smooth initial condition to satisfy the requirement for uq, 
consequently, the theoretical results established in this paper apply to this test prob¬ 
lem. On the other hand, non-smooth initial conditions are used in the second and 
third tests, hence, the theoretical results of this paper may not apply. But we still 
use our MIP-DG methods to compute the error order, energy decay and the evolu¬ 
tion of the numerical interfaces. Our numerical results suggest that the proposed DG 
schemes work well, even a convergence theory is missing for them. 


Testl. Consider the Cahn-Hilliard problem (l.I)-(I.5l with the following initial 
condition: 


uo(x) = tanh 


( 


do(x)\ 

V2e /’ 


where tanh(t) = (e* — e *)/(e* -I- e *), and do(x) represents the signed distance 
function to the ellipse: 


xf 


X2 


0.36 0.04 


= I. 


Hence, uq has the desired form as stated in Proposition 3.5 

Table 5.1 shows the spatial and H^-norm errors and convergence rates, which 


are consistent with what are proved for the linear element in the convergence theorem, 
e = 0.1 is used to generate the table. 



error 

L°°{L^) order 

error 

order 

h = QA^/2 

0.53325 


0.84260 


h = 0.2^2 

0.21280 

1.3253 

0.64843 

0.3779 

h = 0.1^2 

0.07164 

1.5707 

0.43273 

0.5835 

h = 0.05^2 

0.01779 

2.0097 

0.21411 

1.0151 

h = 0.025^2 

0.00454 

1.9703 

0.10890 

0.9753 


Table 5.1 

Spatial errors and convergence rates of Test 1 with e = 0.1. 


Figure 5.1 plots the change of the discrete energy Eh{U^) in time, which should 
decrease according to (3.171. This graph clearly confirms this decay property. Figure 


5.2| displays four snapshots at four fixed time points of the numerical interface with 
four different e. They clearly indicate that at each time point the numerical interface 
converges to the sharp interface Ft of the Hele-Shaw flow as e tends to zero. It also 
shows that the numerical interface evolves faster in time for larger e and confirms the 
mass conservation property of the Cahn-Hilliard problem as the total mass does not 
change in time, which approximates a constant 3.064. 


Test 2. Consider the Cahn-Hilliard problem (1.1)-(1.5| with the following initial 
condition: 

uo{x) = tanh^-^ (min{^(a:i -I- 0.3)^ -I- *2 “ 0.3, f {xi — 0.3)^ -I- — 0.25})^ 
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3.5 
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Fig. 5.1. Decay of the numerical energy of Test 1. 






Fig. 5.2. Test 1: Snapshots of the zero-level set of at time t = 0,0.005,0.015,0.03 and 

e = 0.125,0.025,0.005,0.001. 


We note that uq can be written as 



Here do(ai) represents the signed distance function. We note that uq does not have 
the desire d for m as stated in Proposition |3.5| 

Table shows the spatial and H^-norm errors and convergence rates, which 
are consistent with what are proved for the linear element in the convergence theorem, 
e = 0.1 is used to generate the table. Figure [53| plots the change of the discrete energy 
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error 

L°°{L^) order 

error 

order 

h = 0.4^2 

0.26713 


0.35714 


h = 0.2a/2 

0.07161 

1.8993 

0.18411 

0.9559 

h = 0.1a/2 

0.01833 

1.9660 

0.09620 

0.9365 

h = 0.05a/2 

0.00476 

1.9452 

0.04928 

0.9650 

h = 0.025^2 

0.00121 

1.9760 

0.02497 

0.9808 


Table 5.2 

Spatial errors and convergence rates of Test 2 with e = 0.1. 


Eh{U^) in time, which should decrease according to (3.17). This graph clearly confirms 
this decay property. Figure |5.4| displays four snapshots at four fixed time points of 



Fig. 5.3. Decay of the numerical energy Eh{U^) of Test 2. 


the numerical interface with four different e. They clearly indicate that at each time 
point the numerical interface converges to the sharp interface Tj of the Hele-Shaw 
flow as e tends to zero. It again shows that the numerical interface evolves faster in 
time for larger e and confirms the mass conservation property of the Cahn-Hilliard 
problem as the total mass does not change in time, which approximates a constant 
3.032. 

Tests. Consider the Cahn-Hilliard problem (fn])-(fr^ with the following initial 
condition: 


Uo{x) = tanh^—^ (min{y ^{xi -\- 0.3)^ -I- — 0.2, ^{xi — 0.3)^ 


+ x^-0.2, 


^Jxl + {X2 + 0.3)2 _ 0.2, ^xl + [X2 - 0.3)2 _ 0.2})). 


Notice th at th e above Ug does not have the desired form as stated in Proposition |3.5[ 
Table 5.3 shows the spatial and i7^-norm errors and convergence rates with 
e = 0.1, which are consisten t wi th what are proved for the linear element in the 

plots the change of the discrete energy Eh{U^ 


5.5 


convergence theorem. Figure ^_ 

time, which again decreases as predicted by ( 3.17[ ). Figuredisplays four snapshots 
at four fixed time points of the numerical interface with four different e. Once again, 
we observe that at each time point the numerical interface converges to the sharp 
interface Ft of the Hele-Shaw flow as e tends to zero, the interface evolves faster in 
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Fig. 5.4. Test 2: Snapshots of the zero-level set of at time t = 0,0.001,0.04,0.09 and 

€ = 0.125,0.025,0.005,0.001. 



error 

order 

error 

order 

h = 0.4v^ 

0.38576 


0.84157 


h = 0.272 

0.12347 

1.6435 

0.55082 

0.6115 

h = O.Ia/2 

0.03599 

1.7785 

0.31149 

0.8224 

h = 0.0572 

0.00965 

1.8990 

0.16199 

0.9433 

h = 0.02572 

0.00247 

1.9660 

0.08218 

0.9790 


Table 5.3 

Spatial errors and convergence rates of Test 3 with e = 0.1. 


time for larger e and the mass conservation property is preserved. The total mass 
approximates a constant 2.989. 
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